Tarea 2: Frequentist Inference

CC6104: Statistical Thinking

Integrantes :

Cuerpo Docente:

Fecha límite de entrega:

Índice:

  1. Objetivo
  2. Instrucciones
  3. Referencias
  4. Primera Parte: Preguntas Teóricas
  5. Segunda Parte: Elaboración de Código

Objetivo

a la segunda tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la segunda parte del curso, los cuales se enfocan principalmente en inferencia estadística, diseño de experimentos y test de hipótesis. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

Referencias:

Slides de las clases:

Enlaces a videos de las clases:

Documentación:

Primera Parte: Preguntas Teóricas

A continuación, se presentaran diferentes preguntas que abordan las temáticas vistas en clases. Por favor responda cada una de estas de forma breve.

Pregunta 1:

A continuación, se presenta una serie de declaraciones relacionadas con el sesgo en el muestreo, señale de forma breve el tipo de sesgo que se observa y como podría solucionar el problema de sesgo. Si no observa sesgo en alguna de las declaraciones, comente solamente que es un experimento sin sesgo:

  1. La Tercera quiere conocer la opinión de la gente sobre una propuesta de ley sobre el matrimonio igualitario. Un reportero del diario sale de la central ubicada en Las Condes y selecciona al azar a 300 personas que pasean por allí, preguntándoles sobre la ley propuesta. ¿Qué podemos decir de este plan de muestreo?.

Podemos decir que el plan de muestreo presenta un Sesgo de Selección, debido a que es de común conocimiento que existe mayor densidad de personas con pensamientos políticos más conservadores en la comuna de Las Condes, por lo que es muy probable que gran parte de los encuestados responderá negativamente a esta propuesta de ley. Una forma de solucionar lo anterior sería realizar multiples encuestas en distintas comunas del país y no solamente en Las Condes, dando mayor representatividad al resto de la población.

  1. Una empresa de telefonía quiere conocer el grado de satisfacción de los propietarios de su nuevo servicio 5g, para ver que tal es la calidad del servicio entregado. Para ello, selecciona al azar 500 números entre todo el conjunto de clientes, colocándose en contacto telefónico con ellos. ¿Qué podemos decir de este plan de muestreo?

Este nuevamente es un Sesgo de Selección, debido a que la empresa selecciona al azar entre todo el conjunto de clientes respecto a un producto que seguramente no todos poseen, por lo que no se tendrán respuestas reales respecto al servicio, si no respecto a la imagen que tiene cada cliente. La forma de arreglar este problema es elegir 500 numeros al azar pero del conjunto de clientes que hayan contratado efectivamente el servicio 5g.

  1. Una compañía aérea quiere hacer una encuesta a sus clientes para mejorar su servicio. Durante un mes, envía un correo electrónico a una muestra aleatoria de clientes que hayan volado con la aerolínea el día anterior (ningún cliente será contactado más de una vez). En el correo electrónico se indica que la compañía aérea desea que el cliente rellene una encuesta de 10 minutos para ayudar a la compañía a mejorar su servicio. ¿Qué podemos decir de este plan de muestreo? ¿el sorteo de un producto podría solventar esto?.

Este plan de muestreo tiene un sesgo dificil de reconocer, y es el Sesgo de Respuesta Voluntaria. Debido a que la empresa consulta si el cliente quiere realizar la encuesta, lo más probable es que los clientes altamente satisfechos o clientes altamente descepcionados con el servicio aereo sean propensos a responder la encuesta para plazmar sus opiniones. El sorteo de un producto podría solventar esto, ya que llamaría a más clientes que antes del sorteo no estaban interesados en responder la encuesta.

Pregunta 2:

Explique una buena practica para desarrollar un muestreo, ¿Cree que en una encuesta real es factible obtener una muestra sin sesgo?

Una buena práctica es hacer un Muestreo Aleatorio Estratificado, el cuál consiste en elegir aleatoriamente la muestra pero bajo restricciones de cantidad y proporción respecto a un parámetro específico de la población strata.

Por ejemplo, una encuesta en el Campus Beauchef podría ser Estratificada si se eligiera de este un 40% de mujeres y un 60% de hombres para evitar la sobre representación de hombres en la muestra, donde en este caso la strata es el género. Otra Strata en este mismo ejemplo podría ser la carrera de los alumnos, región de procedencia, etc.

Si hablamos de una encuesta real, y aunque estratifiquemos la muestra, siempre habrá algún tipo de sesgo por muy mínimo que sea y que no se pueda controlar. El mejor ejemplo de esto es el Sesgo de Sampling Frame, debido a que no podemos asegurar con total certeza que el sample frame está hecho a la perfección, o si algún individuo no respondió con total honestidad.

Pregunta 3:

Nombre una ventaja y una desventaja del método de máxima verosimilitud. Justifique cada una de las ventajas y desventajas señaladas.

Una ventaja notable del método de máxima verosimilitud es que es consistente, además de que los estimadores de mazima verosimilitud tienden a una distribución asintótica Gaussiana.

La principal desventaja consiste en que su uso depende de la distribución. En la vida real uno puede adivinar la distribución de situaciones familiares (un lanzamiento de monedas será binomial, la entrada y salida de personas será Poisson, etc), pero para casos más exóticos puede ser difícil determinar esto.

Pregunta 4:

Considere que recopila \(X_{1},...,X_{n} \sim Poisson(\theta)\) donde \(\theta\) es un parámetro desconocido. Calcule el estimador de máxima verosimilitud para \(\theta\). Considere la siguiente expresión para la Poisson:

\[ f(X_{i};\theta) = \prod_{i=1}^{n} \frac{\theta^{X_{i}}e^{-\theta}}{X_{i}!}\]

Nota: Puede ser útil utilizar el logaritmo natural para inferir \(\theta\).

Utilizando verosimilitud logarítmica, tendremos que:

\[ \mathcal{L}_{n}(\theta) = \prod_{i=1}^{n} \frac{\theta^{X_{i}}e^{-\theta}}{X_{i}!} = \frac{\theta^{(\sum_{i=1}^{n} X_{i})}e^{- n \theta}}{\prod_{i=1}^{n} X_{i}!} \]

Y al tomar \(S = \sum_{i=1}^{n} X_{i}\) y aplicar la función logaritmo a la expresión, tenemos que:

\[ l_{n}(\theta) = S \ln(\theta) - n\theta - \sum_{i=1}^{n}\ln(X_{i}!)\]

Donde para inferir el valor de \(\theta\) podemos derivar esta función logarítmica e igual a 0. Así:

\[ \frac{\partial }{\partial \theta} l_{n}(\theta) = \frac{S}{\theta} - n \]

\[\frac{S}{\theta} - n = 0\]

\[ \Rightarrow \theta = \frac{S}{n} = \frac{\sum_{i=1}^{n} X_{i}}{n} = \bar{X}\]

Si \(S\) y \(n\) son distintos de 0.

Pregunta 5:

Suponga que usted ha realizado una sampling distribution de dos encuestas diferentes en los que obtiene los datos Muestreo 1 y Muestreo 2, al momento de estimar el parámetro \(\beta\) de cada uno de los muestreos, se da cuenta que el estimador para los diferentes datos varia según el tamaño de observaciones (n). Al darse cuenta de esto, gráfica el histograma del comportamiento de los estimadores, obteniendo los siguientes gráficos:

Suponiendo que la linea roja en los gráficos representa el valor real del parámetro para cada muestra, comente las características que observa en la evolución del estimador para el Muestreo 1 y el Muestreo 2.

El estimador del Muestreo 1 es claramente consistente, debido a que al ir aumentando el tamaño de la muestra este disminuye su varianza y su media tiende más al valor real.

Por otro lado, pese a que el estimador del Muestreo 2 también disminuye su varianza, tiende a un valor muy distinto al real pero constante en cada iteración. Lo anterior podría indicar que la muestra tiene un sesgo que no depende del tamaño de la muestra, si no a otros parámetros externos como errores metodológicos.

Pregunta 6:

Considere que esta trabajando con intervalos de confianza al \(95\%\) e intenta estudiar un parámetro \(\theta\), ¿porque es incorrecto decir, en la visión frecuentista, que la probabilidad que \(\theta\) pertenezca al intervalo es del \(95\%\)?

Porque un intervalo de confianza no es una afirmación porbabiliística sobre \(\theta\), debido a que el parámetro \(\theta\) es un valor fijo en la visión frecuentista y No una variable aleatoria. La forma correcta de interpretar el intervalo de confianza es que en el \(95\%\) de los experimentos realizados para estimar \(\theta\), el intervalo calculado contiene el valor real de \(\theta\).

Pregunta 7:

Se sabe que al realizar un test de hipótesis es posible equivocarse al tomar una decisión. Sí denotamos \(H_0\) como la hipótesis nula y \(H_1\) como la hipótesis alternativa, existen dos posibles errores que se pueden obtener:

  • Rechazar \(H_0\) cuando \(H_0\) era correcta.
  • Aceptar \(H_0\) cuando \(H_1\) era correcta.

En base a los errores señalados, de un ejemplo de cuando es más relevante disminuir las opciones de aceptar \(H_0\), si \(H_1\) es la correcta.

Cuando estamos estudiando hipótesis relacionadas con la salud de las personas, o dilemas éticos. Por ejemplo:

Queremos realizar un test respecto a si un ascensor resiste menos peso de lo que indica el fabricante. Nuestra hipótesis nula sería que el ascensor resiste lo mismo o más de lo que indica el fabricante, mientras que la hipótesis alternativa es que el ascensor resiste menos de lo que se indica.

En el caso anterior es necesario reducir el margen de error y las posibilidades de aceptar la hipótesis nula cuando la hipótesis alternativa es la correcta, puesto que entra en juego la integridad de las personas que utilizan los ascensores mencionados.

Pregunta 8:

Explique cual es la diferencia entre el test one-sided y two-sided, de un ejemplo en donde se pueda utilizar one-sided y otro donde se pueda utilizar two-sided. ¿Es un método mejor que el otro?

La diferencia entre one-sided y two-sided es que el primer método se utiliza para probar una relación específica con una hipótesis nula (mayor que, menor que), mientras que en el segundo método queremos demostrar una relación de igualdad (mayor o menor, que sea distinto).

Respecto a si un método es mejor que el otro es bastante discutible. El método one-sided es menor restrictivo en cuanto a la hipótesis nula, ya que basta que esté dentro del \(alpha\)% extremo de la hipótesis alternativa para poder demostrarla, mientras que el método two-sided, al desconocer la relación exacta respecto al valor de la hipótesis nula, es doblemente exigente que one-sided al tener que pertenecer al \(alpha/2\)% de cualquier extremo para ser aceptada la hipótesis alternativa.

Pregunta 9:

Explique, con sus palabras, a que se refiere la siguiente afirmación: “Un resultado de significancia estadística es distinto a uno de significancia practica”. Complemente su respuesta analizando la siguiente frase:

'Se tienen datos normales de media desconocida. Con los datos, se realiza un test de hipótesis donde la hipótesis nula especifica que la media es cero. Realizando el test, se obtuvo un p-valor de $10^{-10}$, luego producto de este resultado la media necesariamente tiene que ser muy distinta a cero.'

La afirmación se refiere a que la significancia estadística (como el p-valor, estimador, etc.) no aporta lo mismo en significacia práctica (información útil). En el ejemplo podemos tomar el resultado de significancia estadística como el p-valor de la hipótesis (que en este caso es de \(10^{-10}\)) pero en la práctica no nos sirve saber con exactitud la probabilidad, ya que basta saber que la probabilidad es baja para poder llegar a otras conclusiones relevantes.

Segunda Parte: Elaboración de Código

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 2 son las siguientes:

# Manipulación de estructuras
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.3     v purrr   0.3.4
## v tibble  3.1.0     v dplyr   1.0.5
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   1.4.0     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(Rlab)
## Rlab 2.15.1 attached.
## 
## Attaching package: 'Rlab'
## The following object is masked from 'package:dplyr':
## 
##     count
## The following object is masked from 'package:tibble':
## 
##     view
## The following objects are masked from 'package:stats':
## 
##     dexp, dgamma, dweibull, pexp, pgamma, pweibull, qexp, qgamma,
##     qweibull, rexp, rgamma, rweibull
## The following object is masked from 'package:datasets':
## 
##     precip
# Para realizar plots
library(ggplot2)
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
# Manipulación de varios plots en una imagen.
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.0.5
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library(grid)

Pregunta 1: Estimadores.

Esta pregunta tiene como objetivo comprobar a través de gráficos las características que poseen los estimadores.Por favor responda de forma separada las siguientes preguntas:

  • En clases se vio que el estimador \(\hat{p}_{n} = \frac{1}{n} \displaystyle{\sum_{i=1}^{n}}X_{i}\) es un estimador consistente para \(X_{i}\) Bernoulli de tasa \(p\), verifique esto númericamente para una distribución Bernoulli de \(p=0.5\), es decir grafique como se ve \(\hat{p}_{n}\) para valores de \(n\) y compárelo con el valor verdadero.

  • Sabemos que no todos los estimadores insesgados son consistentes. Considere el estimador \(T_{n} = \hat{p}_{n} + \epsilon_{n}\) donde \(\epsilon_{n} \sim \mathcal{N}(0,1)\) es posible verificar que \(T_{n}\) es insesgado pero no es consistente, para ver esto repita lo que realizo en el punto anterior y estudie lo que sucede. ¿Porque cree que no es consistente el estimador?, Justifique su respuesta.

Respuesta:

Parte 1

Primero, creaemos los datos necesarios para poder generar los gráficos:

# Definición de variables o estructuras necesarias para el muestreo.
prob <- 0.5
m <- c(1:1000)
n_individuos = 5000
means <- c(0)

# Generamos una población que se comporta como Bernoulli(p=0.5)
bernoulli = rbern(n_individuos, p = prob)
real_mean = mean(bernoulli)

# Por cada tamaño de muestras
for(i in m) {
  # Saco una muestra de tamaño i de la población
  sampling <- sample(bernoulli, i, replace = FALSE)
  
  # Las concateno al vector principal y calculo el estimador
  hat_p = mean(sampling)
  means <- append(means, hat_p)
}

Ya obtenidos los datos de los muestreos, pasaremos a graficar las 4 tomas de muestras con sus estimadores (color azul), y además con la visualización del valor real de la media (color rojo):

p_means <- data.frame(values = means)

# Gráfico de la convergencia
graph <- ggplot(p_means, aes(x=c(0:1000), y=values)) + # graficamos el dataframe
  geom_line() +    # graficamos los datos en forma de línea
  ylim(0, 1) +     # ponemos límites a los ejes
  xlim(0, 1000) +  
  geom_hline(aes(yintercept=0.5),  
            color="red", linetype="solid", size=1, alpha = 0.4) +
  xlab("Tamaño Muestra") +  
  ylab("Estimador") + 
  ggtitle("Valores del estimador p en función del tamaño de la muestra") + 
  theme(plot.title = element_text(hjust = 0.5))
graph

Se puede observar una tendencia clara del valor del estimador hasta el valor real de la media a medida que aumentamos el tamaño de cada muestra, indicandonos que este estimador es consistente.

Parte 2

Utilizando la misma idea de la primera parte, crearemos los datos necesarios para pasar a graficarlos posteriormente:

# Definición de variables o estructuras necesarias para el muestreo.
prob <- 0.5
m <- c(1:1000)
n_individuos = 5000
means <- c(0)

# Generamos una población que se comporta como Bernoulli(p=0.5)
bernoulli = rbern(n_individuos, p = prob)
real_mean = mean(bernoulli)

# Por cada tamaño de muestras
for(i in m) {
  # Saco una muestra de tamaño i de la población
  sampling <- sample(bernoulli, i, replace = FALSE)
  
  # Las concateno al vector principal y calculo el estimador
  hat_p = mean(sampling) + rnorm(1, mean=0, sd=1)
  means <- append(means, hat_p)
}

Ya obtenidos los datos de los muestreos, pasaremos a graficar las 4 tomas de muestras con sus estimadores (color azul), y además con la visualización del valor real de la media (color rojo):

p_means <- data.frame(values = means)

# Gráfico de la convergencia
graph <- ggplot(p_means, aes(x=c(0:1000), y=values)) + # graficamos el dataframe
  geom_line() +    # graficamos los datos en forma de línea
  ylim(-2, 2.5) +     # ponemos límites a los ejes
  xlim(0, 1000) +  
  geom_hline(aes(yintercept=0.5),  
            color="red", linetype="solid", size=1, alpha = 0.4) +
  xlab("Tamaño Muestra") +  
  ylab("Estimador") + 
  ggtitle("Valores del estimador p en función del tamaño de la muestra") + 
  theme(plot.title = element_text(hjust = 0.5))
graph

A diferencia del caso anterior, el valor del estimador varía significativamente y de forma aleatoria por cada iteración. Debido a que este no pareciera converger a ningún valor real, y su variabilidad incluye incluso valores fuera de los valores posibles de la distribución, este no es un estimador consistente.

Pregunta 2: Intervalo de Confianza

El objetivo de esta pregunta es visualizar los intervalos de confianza en datos simulados de una población, para visualizar la incertidumbre que presenta una estimación. Para esto, ustedes deberán generar datos de una distribución exponencial, la cual deberán considerar como los datos de la población. En base a los datos generados, estimen la distribución de la media de la población a través de la sampling distribution de la media. Notar que el valor obtenido en cada muestra les entregará un estimador de la media, o sea, para cada valor podremos calcular un intervalo de confianza. Hecho esto, calculen el intervalo de confianza del \(95\%\) para cada una de las medias estimadas, utilizando la función de cuantil vista en clases.

Para la elaboración de esta parte de la tarea, se recomienda realizar el experimento con la siguiente secuencialidad:

  • Obtener la media de la población (en este caso la exponencial).
  • Realizar una sampling distribution con un tamaño de muestra igual a \(30\) sobre los datos generados de la población. Repita la obtención de la media un número elevado de veces (recomendación \(5000\) veces).
  • Calcular el intervalo de confianza del \(95\%\) para cada uno de las medias obtenidas en las iteraciones.
  • De acuerdo a los valores obtenidos (medias e intervalos de confianza), grafique cada una de las medias obtenidas en conjunto a sus intervalos de confianza. Aquí debe notar que, si el intervalo de confianza contiene a la media de la población, este se considerará como parte del intervalo de confianza del \(95\%\), haga un conteo de cuantos valores están contenidos en este intervalo.
  • Señalar en el plot de intervalo de confianza los valores que están dentro y fuera del intervalo. Comente que visualiza de los intervalos de confianza obtenidos, ¿existe incertidumbre?.
  • En base al conteo realizado en el punto anterior, observe cómo se comporta la proporción de intervalos de confianza, ¿es acaso este igual al \(95\%\) teórico usado para calcularlo?.

Notar: Responder cada una de las preguntas señaladas en esta sección.

Hints:

  • Para realizar la sampling distribution podría serle útil el comando sample().
  • Puede ser útil generar un dataframe para graficar con ggplot2.
Gráfico esperado para intervalos de confianza

Del gráfico es posible observar que la línea punteada es la media de la población y los puntos de colores son las estimaciones con sus respectivos intervalos de confianza. Notar que para el plot no se utilizaron las 5000 veces, se recomienda utilizar 100 valores para visualizar bien el fenómeno.

Respuesta:

# Definimos tamaños de muestreo
tamano_muestra = 30
n_muestras = 5000
c_in = 0

# Generamos una exponencial para luego generar el subsampling de ella
exponencial = rexp(10000, rate = 2)
# Obtenemos la media poblacional de la exponencial
exp.mean <- mean(exponencial)

# Creamos el data.frame general para poder graficarlo posteriormente:
confidence <- data.frame(
  mean = double(),
  error = double(),
  group = integer(),
  acum_prob = double() # Lo utilizaremos para plotear las proporciones
)

# Sampling distribution, calculo del intervalo de confianza y proporción.
for (i in 1:n_muestras){
  # En base a la distribución exponencial, generamos multiples sampling distribution.
  sampling <- sample(exponencial, tamano_muestra, replace = FALSE)
  
  # Se estima la media del muestreos y el error estándar
  sample.mean <- mean(sampling)
  sample.sd <- sd(sampling)
  sample.se <- sample.sd/sqrt(tamano_muestra)
  
  # Calculamos el t_score y el margen de error para calcular el intervalo
  alpha = 0.05
  t.score = qt(p=alpha/2, df = tamano_muestra - 1, lower.tail=F)

  margin.error <- t.score * sample.se
  
  # Se determina si la media real pertenece al intervalo o no
  g = "IN"
  
  if (exp.mean >= sample.mean - margin.error && exp.mean <= sample.mean + margin.error){
    # Si la media pertenece al intervalo calculado, es del grupo 0
    g = "IN" 
    c_in = c_in + 1
  } else {
    # Si la media no pertenece al intervalo, es del grupo 1
    g = "OUT"
  }
  
  # Concatenamos al dataframe de la muestra los datos rescatados en cada
  # iteración. Este dataframe contendrá toda la información para plotear.
  z <- data.frame(mean = sample.mean, error = margin.error, group = factor(g), 
                  acum_prob = c_in/i)
  confidence <- rbind(confidence, z)
}
# Plot de Intervalos de confianza (ver respuesta esperada)

# Utilizamos el dataframe creado con cada iteración
plot1 <- ggplot(confidence, aes(mean, c(1:5000), colour = group)) + 
  geom_point() +  # Creamos un gráfico de puntos
  xlim(0, 1.5) +  # Añadimos límites para una visualización comoda
  scale_color_manual(values=c("darkturquoise", "deeppink3")) +  # Colores de cada punto
  # Graficamos con geom_errorbar horizontal el intervalo de confianza por cada
  # iteración. Este se calcula como la media calculada más-menos el error calculado.
  geom_errorbarh(aes(xmax = mean + error, xmin = mean - error)) +
  # Añadimos una recta vertical que indique la media real.
  geom_vline(aes(xintercept=exp.mean),
            color="black", linetype="dotted", size=0.5) +
  # Finalmente añadimos labels y un título.
  labs(x="Means", y="Iteration") +
  ggtitle("Intervalos de Confianza \n n = 5000, t = 30") +
  theme(plot.title = element_text(hjust = 0.5))

plot1

Del gráfico anterior podemos observar una clara diferencia en la cantidad de intervalos clasificados como “IN” de color celeste respecto a los “OUT”, pero la proporción de los datos es incierta debido a la cantidad de datos graficados en este plot, y no es posible observar algún tipo de incertudumbre en particular.

# Plot de proporción de Intervalos de confianza
proportion <- data.frame(values = confidence$acum_prob)

# Gráfico de la convergencia
graph <- ggplot(proportion, aes(x=c(1:5000), y=values)) + # graficamos el dataframe
  geom_line() +    # graficamos los datos en forma de línea
  ylim(0.8, 1) +     # ponemos límites a los ejes
  xlim(1, 5000) +  
  geom_hline(aes(yintercept=0.95),  
            color="red", linetype="solid", size=1, alpha = 0.4) +
  xlab("Iteración") +  
  ylab("Proporción Acumulada") + 
  ggtitle("Proporción de intervalos clasificados como 'IN' por iteración") + 
  theme(plot.title = element_text(hjust = 0.5))
graph

Ya graficada la proporción por cada iteración, es claro que existe algún tipo de incertudumbre que afecta al valor esperado. Podemos observar que al aumentar la cantidad de iteraciones este no converge al valor teórico, si no que a un estimado de 92.5%, indicando que existe algún tipo de incertidumbre que podría depender del tipo de distribución o del tamaño de cada muestra tomada.

Pregunta 3: Estimación de Máxima Verosimilitud

En esta pregunta deberán trabajar con el dataset Body Measurements_original.csv. El objetivo será visualizar e inferir los parámetros que componen a dos variables del dataset. Para esto deberá visualizar el comportamiento de la likelihood, utilizando diferentes cantidades de datos, y realizar la optimización de la likelihood para obtener los estimadores de las variables a través de la función nlminb(). Notar que esta pregunta consiste en dos partes.

  1. La primera actividad consiste en realizar inferencia sobre la variable TotalHeight, donde deberá asumir y realizar los siguientes puntos:
  • Asuma que los datos de TotalHeight distribuyen de forma gaussiana. Visualice esto a través de un gráfico de densidad en la variable TotalHeight.
  • Grafique a través de un gráfico de calor el rango de valores en que se mueve la solución del problema de likelihood, para esto defina su problema en base a -log(likelihood). Mas información se entrega en el esqueleto de la pregunta.
  • Visualice cómo se comporta la -log(likelihood) con diferentes cantidades de datos para la estimación de \(\mu\); para visualizar esto, fije \(\sigma\) en 12 y varié solamente el valor de \(\mu\). Para observar el comportamiento del estimador con diferentes cantidades de datos, realice un subsampling de 100 datos, 300 datos y la totalidad de los datos de la variable TotalHeight. ¿Qué se puede observar acerca del comportamiento del estimador \(\mu\)?, responda brevemente esta pregunta.
  • Finalmente aplique el comando nlminb() sobre la likelihood y encuentre el máximo o mínimo del problema a optimizar. Señale implícitamente cuales son los valores obtenidos para \(\sigma\) y \(\mu\).
  1. Para la segunda parte deberá escoger otra de las variables que componen el dataset, como por ejemplo Age, y estimar a traves de la -log(likelihood) solo los parámetros de la distribución que observa (notar que solo debe inferir los estimadores de variable escogida). Para señalar la distribución de los datos se recomienda realizar el plot de densidad y comparar con el comportamiento de las distribuciones teóricas vistas en clases.

Cabe señalar que el método de máxima verosimilitud deberá ser programado por usted y no podrá utilizar librerías que entreguen el valor directo (por ejemplo, la librería MASS).

Respuesta

# Carga de dataset
body_measurements <- read.csv("Body Measurements _ original.csv")
# Primera Parte
# Plot de densidades de la variables TotalHeight
dens <- density(body_measurements$TotalHeight) # returns the density data
plot(dens) # plots the results

# Plot de Likelihood
normal_pdf <- function(mu,sigma,x) {
  return (1 / (sigma*sqrt(2*pi)) * exp( -1 / 2 * ((x-mu)/sigma)^2 ) )
}
# - Generar una función de la likelihood de la normal

likelihood_normal <- function(mu,sigma,xs) {
  n = length(xs)
  x_barra = mean(xs)
  S2 = sum((xs-x_barra)^2)
           
  return (1/sigma^2)*exp(-(n*S2/(2*sigma^2)))*exp(-(n*(x_barra-mu)^2/(2*sigma^2)))
  #return (2*pi*sigma^2)^(-n/2) * exp( (-1/(2*sigma^2)) * (sum( (xs-mu)^2 )))
  
  
}

# - Señalar el rango de valores para observar la solución. Genere un vector con 
# los valores
valores = as.vector(body_measurements$TotalHeight)

# Plotear gráfico de calor a través filled.contour()
ll_plot <- function(mu, sigma) {
  # Definimos la likelhood
  xs = valores
  return (likelihood_normal(mu,sigma,xs))
}

# Vectorizamos nuestra función para recorrer y estimar los valores
ll_plot = Vectorize(ll_plot)

mu = seq(20,80,by=0.5)# definimos secuencia de 20 -> 80 de 0.5 en 0.5.
sigma = seq(5,23,by=0.5) # definimos secuencia de 5 -> 23 de 0.5 en 0.5.

ll_plot = outer(X=mu, Y=sigma, ll_plot)

# Obtenemos el mapa de calor con los valores mas probables
filled.contour(x=mu, y=sigma, z=ll_plot, xlab=expression(mu), 
               ylab=expression(sigma))

# Plot de comportamiento de SOLO mu al variar la cantidad de datos
logs_100 = c()
logs_300= c()
for (m in mu){
  logs_100 = append(logs_100,-log(likelihood_normal(m,12,sample(valores,100))))
}
for (m in mu){
  logs_300 = append(logs_300,-log(likelihood_normal(m,12,sample(valores,300))))
}

plot(mu,logs_100)
plot(mu,logs_300)


plot(replicate(100, mean(sample(valores, 100, rep=F))),type="l", col="blue", xlab="Iteración de Samples", ylab="mu", main="Variación de mu")
plot(replicate(100, mean(sample(valores, 300, rep=F))),type="l", col="red", xlab="Iteración de Samples", ylab="mu", main="Variación de mu")

# Obtener la solución que minimiza o maximiza la likelihood
# Producto de como funcionan nlminb es necesario definir un nuevo tipo de función
# para encontrar los parametros de la likelihood, por favor revisar estructura entregada.
likelihood <- function(param) {
  # Definimos los parámetros de entrada de la función
  mu = param[1]
  sigma = param[2]
  xs = valores
  n = length(xs) 
  x_barra = mean(xs)
  S2 = sum((xs-x_barra)^2)
           
  # Definimos la likelihood como la suma logaritmica de la función de densidad
  # return (  (-n/2)*log(2*pi) + (-n/2)*log(sigma^2) - (1/(2*sigma^2))*sum((param[3]-mu)^2)  )
  # return (2*pi*sigma^2)^(-n/2) * exp( (-1/(2*sigma^2)) * (sum( (xs-mu)^2 )))
  # return (-log(likelihood_normal(mu,sigma,xs)))
  
  return (-n*log(sigma))-(n*S2/(s*sigma^2))-(n*(x_barra-mu)^2/(1*sigma^2))
  
}
likelihood_100 <- function(param) {
  # Definimos los parámetros de entrada de la función
  mu = param[1]
  sigma = param[2]
  xs = sample(valores,100)
  n = length(xs) 
  x_barra = mean(xs)
  S2 = sum((xs-x_barra)^2)
           
  # Definimos la likelihood como la suma logaritmica de la función de densidad
  # return (  (-n/2)*log(2*pi) + (-n/2)*log(sigma^2) - (1/(2*sigma^2))*sum((param[3]-mu)^2)  )
  # return (2*pi*sigma^2)^(-n/2) * exp( (-1/(2*sigma^2)) * (sum( (xs-mu)^2 )))
  # return (-log(likelihood_normal(mu,sigma,xs)))
  
  return (-n*log(sigma))-(n*S2/(s*sigma^2))-(n*(x_barra-mu)^2/(1*sigma^2))
  
}

likelihood_300 <- function(param) {
  # Definimos los parámetros de entrada de la función
  mu = param[1]
  sigma = param[2]
  xs = sample(valores,300)
  n = length(xs) 
  x_barra = mean(xs)
  S2 = sum((xs-x_barra)^2)
           
  # Definimos la likelihood como la suma logaritmica de la función de densidad
  # return (  (-n/2)*log(2*pi) + (-n/2)*log(sigma^2) - (1/(2*sigma^2))*sum((param[3]-mu)^2)  )
  # return (2*pi*sigma^2)^(-n/2) * exp( (-1/(2*sigma^2)) * (sum( (xs-mu)^2 )))
  # return (-log(likelihood_normal(mu,sigma,xs)))
  
  return (-n*log(sigma))-(n*S2/(s*sigma^2))-(n*(x_barra-mu)^2/(1*sigma^2))
  
}

# Optimizador para encontrar los parametros de la likelihood. Referencia: https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/nlminb
sigma_12= rep(12,length(mu))
nlminb(objective=likelihood, start=c(mu,sigma_12) , lower=c(45,12) , upper=c(52,12) )
nlminb(objective=likelihood_100, start=c(mu,12) , lower=45 , upper=52 )
nlminb(objective=likelihood_300, start=c(mu,12) , lower=45 , upper=52 )

# Segunda Parte
# Graficar la densidad y obtener el parámetro de la variable propuesta.

Pregunta 4: Bootstrap

Para esta pregunta será necesario cargar el dataset SAT_GPA.csv, con el que estudiaremos la correlación entre las variables SAT y GPA. Dentro de las variables: GPA representa el rendimiento académico de un estudiante en el sistema estadounidense, mientras que SAT es una prueba de admisión universitaria en estados unidos. Las actividades por realizar en esta pregunta son las siguientes:

  • Visualizar la correlación de los datos.
  • Obtener la correlación entre los datos de las variables Total, quien representa el resultado de la prueba SAT en USA, y la variable GPA.
  • Programar el método Bootstrap para calcular el error estándar de la correlación. Para la simulación utilice un numero de muestras (B) igual a 5000 o algún otro numero de este orden.
  • Visualizar a través de un gráfico el histograma obtenido al realizar el muestreo con Bootstrap.
  • Obtener el 95% de intervalo de confianza de la estimación por cuantiles como se vio en clases.

Nota: No se permite la utilización de librerías de bootstrap para esta parte.

Respuesta:

sat_gpa <- read.csv("SAT_GPA.csv")
head(sat_gpa)
summary(sat_gpa)
##       Rank           Bangla         English          Botany        Zoology     
##  Min.   :  3.0   Min.   :35.00   Min.   :40.00   Min.   :31.0   Min.   :39.00  
##  1st Qu.:132.2   1st Qu.:56.00   1st Qu.:53.00   1st Qu.:44.0   1st Qu.:51.00  
##  Median :243.5   Median :62.00   Median :57.00   Median :50.5   Median :55.00  
##  Mean   :244.0   Mean   :62.72   Mean   :56.75   Mean   :49.4   Mean   :54.37  
##  3rd Qu.:360.8   3rd Qu.:70.00   3rd Qu.:61.00   3rd Qu.:55.0   3rd Qu.:59.00  
##  Max.   :520.0   Max.   :84.00   Max.   :72.00   Max.   :65.0   Max.   :65.00  
##                  NA's   :1                                                     
##     Physics        Chemistry          Math           Ict            Total      
##  Min.   :31.00   Min.   :30.00   Min.   :13.0   Min.   :46.00   Min.   :305.0  
##  1st Qu.:45.00   1st Qu.:45.00   1st Qu.:29.0   1st Qu.:56.00   1st Qu.:398.0  
##  Median :51.00   Median :50.00   Median :33.5   Median :60.00   Median :420.5  
##  Mean   :49.79   Mean   :50.40   Mean   :34.2   Mean   :59.53   Mean   :417.2  
##  3rd Qu.:56.00   3rd Qu.:55.75   3rd Qu.:39.0   3rd Qu.:64.00   3rd Qu.:441.0  
##  Max.   :65.00   Max.   :71.00   Max.   :58.0   Max.   :72.00   Max.   :492.0  
##                                                                                
##       Gpa       
##  Min.   :2.500  
##  1st Qu.:3.580  
##  Median :3.920  
##  Mean   :3.897  
##  3rd Qu.:4.200  
##  Max.   :5.000  
## 
sapply(sat_gpa, function(x) sum(is.na(x)))
##      Rank    Bangla   English    Botany   Zoology   Physics Chemistry      Math 
##         0         1         0         0         0         0         0         0 
##       Ict     Total       Gpa 
##         0         0         0
sat_gpa <- sat_gpa[!is.na(sat_gpa$Bangla),]
sapply(sat_gpa, function(x) sum(is.na(x)))
##      Rank    Bangla   English    Botany   Zoology   Physics Chemistry      Math 
##         0         0         0         0         0         0         0         0 
##       Ict     Total       Gpa 
##         0         0         0
# Calculo de la matriz de correlacion para su visualización
pearson_matrix = cor(sat_gpa, method = "pearson") 
round(pearson_matrix, 3)
##             Rank Bangla English Botany Zoology Physics Chemistry   Math    Ict
## Rank       1.000 -0.286  -0.383 -0.667  -0.522  -0.620    -0.602 -0.514 -0.514
## Bangla    -0.286  1.000   0.051  0.256   0.056   0.072    -0.164 -0.038  0.150
## English   -0.383  0.051   1.000  0.203   0.252   0.101     0.256  0.071  0.283
## Botany    -0.667  0.256   0.203  1.000   0.498   0.342     0.264  0.169  0.326
## Zoology   -0.522  0.056   0.252  0.498   1.000   0.226     0.274  0.141  0.123
## Physics   -0.620  0.072   0.101  0.342   0.226   1.000     0.421  0.303  0.214
## Chemistry -0.602 -0.164   0.256  0.264   0.274   0.421     1.000  0.274  0.238
## Math      -0.514 -0.038   0.071  0.169   0.141   0.303     0.274  1.000  0.157
## Ict       -0.514  0.150   0.283  0.326   0.123   0.214     0.238  0.157  1.000
## Total     -0.945  0.383   0.464  0.702   0.546   0.624     0.581  0.497  0.544
## Gpa       -0.981  0.309   0.400  0.672   0.525   0.622     0.595  0.552  0.511
##            Total    Gpa
## Rank      -0.945 -0.981
## Bangla     0.383  0.309
## English    0.464  0.400
## Botany     0.702  0.672
## Zoology    0.546  0.525
## Physics    0.624  0.622
## Chemistry  0.581  0.595
## Math       0.497  0.552
## Ict        0.544  0.511
## Total      1.000  0.966
## Gpa        0.966  1.000
# Visualización de la matriz de correlación
library(corrplot)
## corrplot 0.90 loaded
corrplot(pearson_matrix, type = "upper", tl.col = "black", tl.srt = 50)

# Valor de la correlación de Total y GPA:
pearson_matrix[11, 10]
## [1] 0.9655365
# Utilizado la función de ejemplo de bootstrap visto en clases
mybootstrap <- function(x, fun, nRuns, size, alpha){
  values <- vector()
  for (i in 1: nRuns){
    samp.i <- sample(x, size = size, replace = T)
    values[i] <- fun(samp.i)
  }
  
  point.est <- fun(x)
  se <- sd(values)
  l.CI <- quantile(values, alpha/2)
  u.CI <- quantile(values, 1-alpha/2)
  
  return(c("Point Estimate" = point.est,
           "Standart error" = se,
           "Lower CI Limit" = l.CI,
           "Upper CI Limit" = u.CI))
}
mybootstrap(sat_gpa[10:11], cor , 5000, 125, 0.05)
# Utilizado la función de ejemplo de bootstrap visto en clases
mybootstrap <- function(params, fun, nRuns, size, alpha){
  values <- vector()
  samp <- vector()
  for (i in 1: nRuns){
    samp.i <- sample_n(params, size=size, replace = T)
    values[i] <- fun(samp.i)[2]
  }
  hist(values,breaks="FD")
  point.est <- fun(params)
  se <- sd(values)
  l.CI <- quantile(values, alpha/2)
  u.CI <- quantile(values, 1-alpha/2)
  
  return(c("Point Estimate" = point.est,
           "Standart error" = se,
           "Lower CI Limit" = l.CI,
           "Upper CI Limit" = u.CI))
}
mybootstrap(sat_gpa[10:11], cor , 5000, 125, 0.05)

##      Point Estimate1      Point Estimate2      Point Estimate3 
##          1.000000000          0.965536466          0.965536466 
##      Point Estimate4       Standart error  Lower CI Limit.2.5% 
##          1.000000000          0.006086894          0.952009159 
## Upper CI Limit.97.5% 
##          0.975949193
 

A work by CC6104